Introduction¶
In this case study, we will be working with missing data using the Boston house-price data. We will simply be exploring the change in results by performing imputations for the missing values for multiple variables.
import numpy as np
import pandas as pd
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import accuracy_score
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt
from pandas_profiling import ProfileReport
from IPython.display import Markdown as md
import seaborn as sns
def take_two(data,
null_cols=[],
test_col='AGE',
cond='>10',
impute_func=np.mean,
samples=[.1, .2, .3]):
assert callable(impute_func)
assert isinstance(test_col, str)
assert isinstance(cond, str)
data_frames = {}
for sample_rate in samples:
df = data.copy()
replace_index = df.query(f'{test_col}{cond}').sample(
frac=sample_rate, replace=False).index
print(f'Percent of imputed values:{(len(replace_index)/len(df.index)):,.2%}')
df.loc[replace_index, null_cols] = None
imputed = impute_func(df.loc[df.index.difference(replace_index),
null_cols],
axis=0)
print(f'{imputed}\n{impute_func}\n{sample_rate}', end='\n')
[
df[x].fillna(imputed[i], inplace=True)
for i, x in enumerate(null_cols)
]
data_frames.update({sample_rate: df})
return data_frames
# load data
boston = load_boston()
# display data descriptions
md(boston.DESCR)
DataFrame Creation & Data Partitioning¶
# create independent and dependent dataframes
df = pd.DataFrame(boston.data, columns=boston.feature_names) # features dataframe
df_targets = pd.DataFrame(boston.target, columns=['target']) # targets dataframe
# partition data
X_train, X_test, y_train, y_test = train_test_split(df, df_targets.target, test_size=0.33, random_state=7337)
Explore Data¶
Feature Correlation¶
fig, ax = plt.subplots(figsize=(12,12))
im = ax.matshow(df.corr())
ax.set_xticks(np.arange(len(df.columns))) # force to show all x labels
ax.set_yticks(np.arange(len(df.columns))) # force to show all y labels
ax.set_xticklabels(df.columns)
ax.set_yticklabels(df.columns)
fig.colorbar(im)
plt.show()
Based on the correlation matrix above, we chose to interpret the following variables:
B - 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town LSTAT - % lower status of the population Age - proportion of owner-occupied units built prior to 1940
- LSTAT has a positive correlation with CRIM, as the crime rate by town(CRIM) increases relates to the lower stataus of the population in that part of town There is also a positive correlation between LSTAT and AGE. The percent of lower status of the population is high in areas where proportion of owner-occupied units built prior to 1940 are also high. Typically, neighborhoods with old housing tend to have a lower status of population since they are typically cheaper. Where as newer homes tend to cater to those with higher status of the population because they are able to afford this type of housing.
- Conversely, B and AGE are inversely correlated since proportion of homes built prior to 1940 tend to have a low percentage of proportion of blacks by town. For example, neighborhoods that were prestigious prior to 1940 and still carry that prestige today still have a low percentage of proportion of black people living in these neighborhoods.
- LSTAT is inversely correlated to ZN, RM, DIS. As the % of lower status of the population increases, the ZN, RM, DIS drop. For example, the average number of rooms in houses(RM) in lower status neighborhoods is low compared to number of rooms in houses in more prestigious neighborhoods.
Profile Dataset¶
profile = ProfileReport(df, title="Boston Features Profile")
profile.to_notebook_iframe()
Step 2:¶
Randomize 1%, 5% 10%, 20%, 33%, and 50% missing values in a single column
df_dict = {} # will hold dataframes with various % of missing values
percentages = [0.01, 0.05, 0.1, 0.2, 0.33, 0.5]
column = 'AGE'
for p in percentages:
in_sample = df.sample(frac=p, random_state=42) # sample copy dataframe with fraction at p
in_sample[column] = np.NaN # set selected column values to NaN
out_sample = df[~df.isin(in_sample)].dropna()
df_dict[str(int(p * 100))] = pd.concat([out_sample, in_sample]).sort_index()
for key in df_dict.keys():
print(f"df_dict[\"{key}\"]: '{column}' column, in {int(key) / 100}%",
f"sampled dataframe, has {df_dict[key].isnull().sum()[column]}",
f"of {df_dict[key].shape[0]} missing values")
Baseline Model w/o Missing Values¶
print("Baseline (zero NAs):")
reg = LinearRegression(normalize=True, fit_intercept=True)
reg.fit(X_train, y_train)
#R^2
print("R^2 (Train):", reg.score(X_train, y_train))
print("R^2 (Test):", reg.score(X_test, y_test))
# MSE
pred = reg.predict(X_test)
print("MSE:", mean_squared_error(y_test, pred))
Let's get the fitted values from the model, since sk-learn sucks and doesn't give them to you in the call to .fit()
fitted = fit.predict(X_train)
Plot of Fitted Values v Actual (baseline)¶
fig, ax = plt.subplots(figsize=(20, 10), nrows=1, ncols=1)
p = pd.DataFrame({
'index': range(len(fitted)),
'predicted_values': fitted,
'actual_values': y_train
}).set_index('index').plot(kind='line', ax=ax)
p.lines[0].set_linestyle('--')
p.set_xlabel('Observation')
p.set_ylabel('Median Home Price')
p.set_title('Baseline Fitted vs Actual')
Question 2: Impute data and compare¶
Compare loss and goodness of fit to your baseline
print('Imputation with median value')
print()
for key in df_dict.keys():
print(key, '%')
print(f"Synthetic missing values (at random) for {key}% of '{column}' column.")
_X_train, _X_test, _y_train, _y_test = train_test_split(df_dict[key].copy(), df_targets.target, test_size=0.33, random_state=7337)
# calculate median AGE for imputation
median = _X_train[column].median()
print(f"Imputed with median, {median}")
# impute NAs
_X_train[column].fillna(median, inplace=True)
# build and fit model
reg = LinearRegression(normalize=True)
reg.fit(_X_train, _y_train)
# imput validation set (test) with same value
_X_test[column].fillna(median, inplace=True)
# R^2
print("R^2:", reg.score(_X_test, _y_test))
# MSE
pred = reg.predict(_X_test)
print("MSE:", mean_squared_error(_y_test, pred))
print('-'*20)
print()
Compared to the baseline, the six models with varying degree of missing values from the 'AGE' column had slightly worse MSE and $R^2$ when imputing the newly introduced NAs with the median value for that column. The median value was selected over the mean due to its robustness against outlier values. That being said, mean was also tested with little effect.
Imputation seemed to be a successful approach to achieve comparable losses to baseline model with no missing values in 'AGE'.
Baseline Parameters
Unfortunately sklearn is not R and does not give confidence intervals for parameter estimates, but we can still see the effect that imputation has on the structure of the model by comparing how the parameter estimates change as a result of our imputation scheme.
keys = list(X_train.columns)
keys.append('Intercept')
params = list(fit.coef_)
params.append(fit.intercept_)
baseline_params = pd.DataFrame(zip(keys, params)).rename(columns={
0: 'Field',
1: 'Baseline Parameter Estimate'
})
baseline_params
Step 3:¶
Take two columns and create data “Missing at Random” when controlled for a third variable (i.e., if Variable Z is > 30, then Variables X, Y are randomly missing). Use your preferred imputation method to fill in 10%, 20% and 30% of your missing data.
The function defined below will take any arbitrary combination of columns to be imputed based on a criterion specified for another column. An impution function will be applied to a set of records based on the number of records that meet the criterion and a sample rate. This could be used for a much more exhaustive search for impution schemes if one were so inclined.
For this test - B and LSTAT will be imputed to the mean wherever AGE is less than 40.
samples = take_two(X_train,
null_cols=['B', 'LSTAT'],
test_col='AGE',
cond='<60',
samples=np.linspace(0, 1, 100),
impute_func=np.mean)
X_train[X_train.AGE < 60].shape[0]/X_train.shape[0]
Setting a 10/20/30% sample rate for the rows that match our criteria equates to imputing 3.83%/7.67%/11.21% of our values to the target scheme (in this case, the mean). Note that the mean for B declines as the sample rate increases (354.18 - 351.3) whereas LSTAT increases as the sample rate increases (12.84-13.26).
Refitting LR model to new data frames¶
for rate, data in samples.items():
print(f'{rate:.2%} sample rate results:')
reg = LinearRegression(normalize=True, fit_intercept=True)
fit = reg.fit(data, y_train)
train_score = reg.score(X_train, y_train)
test_score = reg.score(X_test, y_test)
print(
f'Train: {train_score:.6f} is {" better than " if train_score > baseline_score else " worse than " }baseline R2'
)
print(
f'Test: {test_score:.6f} is {" better than " if test_score > baseline_score_test else " worse than " }baseline R2'
)
Question 3:¶
In each case [10%, 20%, 30%] perform a fit with the imputed data and compare the loss and goodness of fit to your baseline. [Note: you should have (9) models to compare against your baseline at this point.]
Using the mean impution scheme described above, all 3 models were slightly better than our baseline model on the testing set - on the training set, they were all slightly worse. The best of the three models (out of sample) was the 30% sample rate model (0.75696 vs 0.75226 $R^2$)
For fun, let's try the same example with the median
samples = take_two(X_train,
null_cols=['B', 'LSTAT'],
test_col='AGE',
cond='<60',
impute_func=np.median)
for rate, data in samples.items():
print(f'{rate:.2%} sample rate results:')
reg = LinearRegression(normalize=True, fit_intercept=True)
fit = reg.fit(data, y_train)
train_score = reg.score(X_train, y_train)
test_score = reg.score(X_test, y_test)
print(
f'Train: {train_score:.6f} is {" better than " if train_score > baseline_score else " worse than " }baseline R2'
)
print(
f'Test: {test_score:.6f} is {" better than " if test_score > baseline_score_test else " worse than " }baseline R2'
)
For fun, let's try the same example with the median but also change the cutoff from 60 years to 85 years. We wrote the code to be experimental after all!
samples = take_two(X_train,
null_cols=['B', 'LSTAT'],
test_col='AGE',
cond='<85',
impute_func=np.median)
for rate, data in samples.items():
print(f'{rate:.2%} sample rate results:')
reg = LinearRegression(normalize=True, fit_intercept=True)
fit = reg.fit(data, y_train)
train_score = reg.score(X_train, y_train)
test_score = reg.score(X_test, y_test)
print(
f'Train: {train_score:.6f} is {" better than " if train_score > baseline_score else " worse than " }baseline R2'
)
print(
f'Test: {test_score:.6f} is {" better than " if test_score > baseline_score_test else " worse than " }baseline R2'
)
Using the median impution scheme described above, all 3 models were slightly better than our baseline model on the testing set - on the training set, they were all slightly worse. The best of the three models (out of sample) was the 30% sample rate model (0.759534 vs 0.75226 $R^2$). This particular scheme (median/>85 cutoff) produced the best results in this analysis. Note that this also had the most samples imputed - 17.40% of the data.
Step 4:¶
Create a “Missing Not at Random” pattern in which 25% of the data is missing for a single column.
Let's cut off the top 25% of AGE since it seems to be very skewed and then impute the missing data to mean.
cuts, bins = pd.qcut(X_train['AGE'], q=[0, .25, .5, .75, 1.], labels=['Q1', 'Q2', 'Q3', 'Q4'], retbins=True)
bins
Cutting out the top 25% of the Age range leaves houses that are less than 94.05 years old.
sample = X_train.copy()
remove_index = cuts[cuts=='Q4'].index
sample.loc[remove_index, 'AGE'] = None
sample_mean = sample.loc[sample.index.difference(remove_index), 'AGE'].mean()
sample.loc[remove_index, 'AGE'] = sample_mean
sample['AGE'].describe()
Question 4:¶
Perform a fit with the imputed data [25%] and compare the loss and goodness of fit to your baseline. [Note: you should have (10) models to compare against your baseline at this point.]
reg = LinearRegression(normalize=True, fit_intercept=True)
fit = reg.fit(sample, y_train)
train_score = reg.score(X_train, y_train)
test_score = reg.score(X_test, y_test)
print(
f'Train: {train_score:.6f} is {" better than " if train_score > baseline_score else " worse than " }baseline R2'
)
print(
f'Test: {test_score:.6f} is {" better than " if test_score > baseline_score_test else " worse than " }baseline R2'
)
Step 5:¶
Describe your imputation approach and summarize your findings. What impact did the missing data have on your baseline model’s performance?
Baseline Score = 0.75226
In Step 3, we decided to impute variables "B" and "LSTAT" to the mean whererver "AGE" is less than 60. The take_two function will take any arbitrary combination of columns to be imputed based on a criterion specified for another column. The specified impution function will be applied to a set of records based on the number of records that meet the criterion and a random sample based on the sample rate. Using the mean imputation scheme described above, all 3 models (in step 3) were slightly better than our baseline model on the testing set - on the training set, they were all slightly worse. The best of the three models (out of sample) was the 30% sample rate model (0.75696 vs 0.75226 $𝑅^2$ )
Given the definitions of these variables (chosen specifically for their social context - see data dictionary) the simplest interpretation of these trends is that as we impute more values of homes with ages less than 60, the concetration of African Americans decreases - this imples that older homes (>60) are located in areas with lower concentrations of African Americans. Conversely, the concentration of those considered "lower-class" increased as we imputed values of homes with ages less than 60. Again, this points to a higher concetration of lower-class individuals with homes over 60 years in age.
After experimenting with the mean imputation scheme, we decided to try the same scheme using the median[10%, 20%, 30%]. We found that all 3 models were slightly better than our baseline model on the testing set - on the training set, they were all slightly worse. The best of the three models (out of sample) was the 30% sample rate model (0.759534 vs 0.75226 $𝑅^2$ ). This particular scheme (median/>85 cutoff) produced the best results in this analysis. Note that this also had the most samples imputed - 17.40% of the data.
In step 4, we decided to cut off the top 25% of AGE since it seems to be very skewed and then impute the missing data to mean. Cutting out the top 25% of the Age range leaves houses that are less than 94.05 years old.